library(coin)
library(dplyr)
library(mosaic)
library(foreign)
library(ggplot2)

filepath <- ('YOUR_FILE_PATH/respondent_game_survey.csv')

data <- read.csv(filepath)[-1]


#dummy out treatment
data$ctr <- ifelse(data$treatment==3, 1, 0)
data$trt2 <- ifelse(data$treatment==2, 1, 0)
data$trt1 <- ifelse(data$treatment==1, 1, 0)


ctr <- data[which(data$treatment==3), ]
trt1 <- data[which(data$treatment==1), ]
trt2 <- data[which(data$treatment==2), ]


# Difference in means + st. dev
means <- c(mean(ctr$payout),mean(trt1$payout,na.rm=T),mean(trt2$payout,na.rm=T))
sqrt(c(var(ctr$payout),var(trt1$payout,na.rm=T),var(trt2$payout,na.rm=T)))

with(data, t.test(payout[data$trt1==1], payout[data$trt1==0]),alternative="greater")
with(data, t.test(payout[data$trt2==1], payout[data$trt2==0]),alternative="greater")

dens <- data[which(complete.cases(data$payout==T)), ]
plot(density(dens$payout),main="Overall and group densities",xlab="Payout")
lines(density(dens$payout[which(dens$trt1==1)]),col="red")
lines(density(dens$payout[which(dens$trt2==1)]),col="blue")
legend("topright",legend=c("overall","horizontal","vertical"),col=c("black","red","blue"),lty=c(1,1))





# Permutation test
numsim = 1000

perm.sub <- c("treatment", "payout", "trt1", "trt2") 
perm.dat <- data[perm.sub]

par(mfrow=c(2,1))

# Treatment 1
res1 = do(numsim) * lm(payout ~ shuffle(trt1), data=data)
obsdiff1 = means[2] - means[1]
pvalue1 = sum(abs(res1$trt1) > abs(obsdiff1)) / numsim

s <- sample(perm.dat$trt1, length(perm.dat$trt1), FALSE)
diff(by(perm.dat$payout, s, mean))

dist <- replicate(2000, diff(by(perm.dat$payout, 
                                sample(perm.dat$trt1, 
                                       length(perm.dat$trt1), 
                                       FALSE), 
                                mean)))

hist(dist, col = "grey", breaks = 100, main="Horizontal",xlab="Difference in Means")
abline(v = diff(by(perm.dat$payout, perm.dat$trt1, mean)), col = "blue", lwd = 2)

# Treatment 2
res2 = do(numsim) * lm(payout ~ shuffle(trt2), data=data)
obsdiff2 = means[3] - means[1]
pvalue2 = sum(abs(res2$trt2) > abs(obsdiff2)) / numsim

s <- sample(perm.dat$trt2, length(perm.dat$trt2), FALSE)
diff(by(perm.dat$payout, s, mean))

dist <- replicate(2000, diff(by(perm.dat$payout, 
                                sample(perm.dat$trt2, 
                                       length(perm.dat$trt2), 
                                       FALSE), 
                                mean)))

hist(dist, col = "grey", breaks = 100,main="Vertical", xlab="Difference in Means")
abline(v = diff(by(perm.dat$payout, perm.dat$trt2, mean)), col = "blue", lwd = 2)


